/* =========================================================== */
/* Module gaseous_transmittance.c                              */
/*                                                             */
/* Computes sensor-specific transmittance through various      */
/* atmospheric gases.                                          */
/*                                                             */
/* B. Franz, NASA/OBPG, July 2006                              */
/* =========================================================== */

#include "l12_proto.h"

void ozone_transmittance(l1str *l1rec, int32_t ip) {
    float tau_oz;
    int32_t iw;

    filehandle *l1file = l1rec->l1file;
    int32_t nwave = l1file->nbands;
    int32_t ipb = ip*nwave;

    for (iw = 0; iw < nwave; iw++) {
        tau_oz = l1rec->oz[ip] * l1file->k_oz[iw];
        l1rec->tg_sol[ipb + iw] *= exp(-(tau_oz / l1rec->csolz[ip]));
        l1rec->tg_sen[ipb + iw] *= exp(-(tau_oz / l1rec->csenz[ip]));
    }
}

void co2_transmittance(l1str *l1rec, int32_t ip) {
    static float *t_co2 = NULL;
    int32_t iw;

    filehandle *l1file = l1rec->l1file;
    int32_t nwave = l1file->nbands;
    int32_t ipb = ip*nwave;

    if (t_co2 == NULL) {
        rdsensorinfo(l1file->sensorID, input->evalmask, "t_co2", (void **) &t_co2);
    }

    for (iw = 0; iw < nwave; iw++) {
        l1rec->tg_sol[ipb + iw] *= pow(t_co2[iw], 1.0 / l1rec->csolz[ip]);
        l1rec->tg_sen[ipb + iw] *= pow(t_co2[iw], 1.0 / l1rec->csenz[ip]);
    }
}

void no2_transmittance(l1str *l1rec, int32_t ip) {

    float a_285, a_225;
    float tau_to200;
    float no2_tr200;
    int32_t iw;

    filehandle *l1file = l1rec->l1file;
    int32_t nwave = l1file->nbands;
    int32_t ipb = ip*nwave;

    float sec0 = 1.0 / l1rec->csolz[ip];
    float sec = 1.0 / l1rec->csenz[ip];

    if (l1rec->no2_tropo[ip] > 0.0)
        /* compute tropo no2 above 200m (Z.Ahmad)    
        no2_tr200 = exp(12.6615 + 0.61676*log(no2_tropo));
           new, location-dependent method */
        no2_tr200 = l1rec->no2_frac[ip] * l1rec->no2_tropo[ip];
    else
        no2_tr200 = 0.0;


    for (iw = 0; iw < nwave; iw++) {

        if (l1file->k_no2[iw] > 0.0) {

            a_285 = l1file->k_no2[iw] * (1.0 - 0.003 * (285.0 - 294.0));
            a_225 = l1file->k_no2[iw] * (1.0 - 0.003 * (225.0 - 294.0));

            tau_to200 = a_285 * no2_tr200 + a_225 * l1rec->no2_strat[ip];

            l1rec->tg_sol[ipb + iw] *= exp(-(tau_to200 * sec0));
            l1rec->tg_sen[ipb + iw] *= exp(-(tau_to200 * sec));

        }
    }
}

void h2o_transmittance(l1str *l1rec, int32_t ip) {
    static float *a_h2o = NULL;
    static float *b_h2o = NULL;
    static float *c_h2o = NULL;
    static float *d_h2o = NULL;
    static float *e_h2o = NULL;
    static float *f_h2o = NULL;
    static float *g_h2o = NULL;

    float t_h2o;
    int32_t iw;

    filehandle *l1file = l1rec->l1file;
    int32_t nwave = l1file->nbands;
    int32_t ipb = ip*nwave;
    float wv = l1rec->wv[ip];

    if (a_h2o == NULL) {
        rdsensorinfo(l1file->sensorID, input->evalmask, "a_h2o", (void **) &a_h2o);
        rdsensorinfo(l1file->sensorID, input->evalmask, "b_h2o", (void **) &b_h2o);
        rdsensorinfo(l1file->sensorID, input->evalmask, "c_h2o", (void **) &c_h2o);
        rdsensorinfo(l1file->sensorID, input->evalmask, "d_h2o", (void **) &d_h2o);
        rdsensorinfo(l1file->sensorID, input->evalmask, "e_h2o", (void **) &e_h2o);
        rdsensorinfo(l1file->sensorID, input->evalmask, "f_h2o", (void **) &f_h2o);
        rdsensorinfo(l1file->sensorID, input->evalmask, "g_h2o", (void **) &g_h2o);
    }

    for (iw = 0; iw < nwave; iw++) {
        t_h2o = a_h2o[iw] + wv * (b_h2o[iw] + wv * (c_h2o[iw] + wv * (d_h2o[iw]
                + wv * (e_h2o[iw] + wv * (f_h2o[iw] + wv * g_h2o[iw])))));
        l1rec->tg_sol[ipb + iw] *= pow(t_h2o, 1.0 / l1rec->csolz[ip]);
        l1rec->tg_sen[ipb + iw] *= pow(t_h2o, 1.0 / l1rec->csenz[ip]);
    }
}


//(int gasmask, int32_t sensorID, int32_t evalmask, int32_t nwave, float mu0, float mu, 
//               float ozone, float wv, float no2_tropo, float no2_strat, float no2_frac,
//	       float t_gas_sol[], float t_gas_sen[])

void gaseous_transmittance(l1str *l1rec, int32_t ip) {
    int32_t iw;

    filehandle *l1file = l1rec->l1file;
    int32_t ipb = ip * l1file->nbands;


    if ((input->gas_opt & O3_BIT) != 0) {
        ozone_transmittance(l1rec, ip);
    }

    if ((input->gas_opt & CO2_BIT) != 0) {
        co2_transmittance(l1rec, ip);
    }

    if ((input->gas_opt & NO2_BIT) != 0) {
        no2_transmittance(l1rec, ip);
    }

    if ((input->gas_opt & H2O_BIT) != 0) {
        h2o_transmittance(l1rec, ip);
    }
}
